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I. INTRODUCTION 

The anchoring phenomenon is the tendency of a liquid 
crystal to orient in a particular direction when in contact 
with the container walls. The equilibrium director orien- 
tation set by the interaction of the liquid crystal with the 
aligning surface is called an easy orientation axis. The 
simplest, strong anchoring, assumption is that the direc- 
tor has a fixed orientation at the boundaries along the 
easy orientation axis. However, it has been discovered 
that the coupling of the director with the orienting sur- 
face can be rather weak. This results in deviation of the 
surface director from the easy axis in response to small 
perturbations. 

On a phenomenological level, weak anchoring can be 
described by adding an appropriate surface potential to 
the free energy of the system. Then minimization of the 
free energy functional gives both the equations for the 
director in the bulk and the appropriate boundary con- 
ditions [0. The simplest form of the surface potential 
has been proposed by Rapini and Papoular B 



f s = - l -W(n-e)\ 



(1) 



where the parameter, W, is termed an anchoring energy. 

Since then, numerous experimental methods have been 
used to measure the surface anchoring coefficient W 
j| ||, ||; its value has turned out to be extremely im- 
portant for liquid crystal devices, i.e. displays, optical 
switches. However, comparatively little work on sys- 
tematic experimental investigation of the anchoring phe- 
nomenon has been presented up till now. From the avail- 
able experimental data, one can say that the extrapola- 
tion length A = is inversely proportional to the 
squared value of the bulk order parameter A cx Q~ 2 . Tak- 
ing into account that the elastic constant -K33 is typically 
proportional to Q 2 , this gives for the anchoring parame- 
ter W oc Q 4 §. 

There have been several attempts to estimate the an- 
choring coefficient theoretically Q and by combining 
molecular simulation with a local density functional ap- 
proach g |. The main difficulty here is that one needs 



to know the direct pair correlation function of the ne- 
matic state, which is usually unknown, and must hence 
be estimated with some uncontrolled approximations or 
assumptions (see, however, Ref. [0). Moreover, it is of- 
ten assumed that the director in the interface region is 
fixed ||. However, it has already been noticed that these 
approximations may give incorrect results, by at least an 
order of magnitude, for example in the calculation of bulk 
elastic constants 
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Similar studies have also been done using lattice mod- 
els. The existence of subsurface deformations and ef- 
fective intrinsic anchoring arising from the incomplete 
molecular interaction close to the surface has been stud- 
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ied using a hexagonal lattice approach 11 1 . Monte-Carlo 
simulations of the Lebwohl-Lasher model have shown 
that the extrapolation length is not in general equal to 
the ratio of the bulk to surface couplings fl^| . However, 
the results were not sufficiently robust to determine the 
dependence of the anchoring energy on the order param- 
eter. 

The present work attempts to remedy the situation 
by combining the elastic description with Onsager the- 
ory and Monte Carlo simulation results. We study the 
dependence of the elastic (K 33 ) and surface anchoring 
(W) coefficients on the liquid crystal state point, which 
is defined by the bulk value of the order parameter Q. 

The paper is organized as follows. We define the ge- 
ometry and derive direc tor profiles using elastic theory 
in section [n| In section III we discuss the Onsager ap- 
proach which allows us to calculate the single particle 
density of the bulk and confined systems, while section 
IV outlines the fluctuation approach to calculating elas- 
tic coefficients. Section |V| gives details of the technique 
used in Monte Carlo simulations and the results are sum- 
marized in section |Vl]. 

II. ELASTIC DESCRIPTION 

The easiest way to obtain the surface anchoring coeffi- 
cient W is to create a director deformation far from the 
surface. Then, measuring the response of the director 
near the cell surface, or fitting the director profile with a 
theoretically predicted profile, yields an estimate of the 
anchoring extrapolation length and ratios of the elastic 
constants. 
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Indeed, consider one of the possible geometries suitable 
for the measurement of the zenithal anchoring strength. 
Let the director have fixed orientation at the boundary 
z = L. The surface at z = is assumed to provide 
homeotropic anchoring of strength W . In the elastic 
description, deformations of the director field n are de- 
scribed by the total free energy pi 



Ta = fbdV+ f s dS. 
Iv Js 



(2) 



Here /b is the Frank-Oseen elastic free energy density 

fb = \ {Kix (V • nf + K 22 (n • [V x n]) 2 + (3) 
+K 33 [nx [Vxn]] 2 }, 

where Kn, K22, and K 33 are elastic constants, and the 
integration extends over the sample volume V. f s is the 
surface anchoring energy density, which we assume to be 
of the Rapini-Papoular form (Q), and it is integrated over 
the boundary surface S. 

In slab geometry, the director n is assumed to lie in the 
xz plane and depends only on the z coordinate. Then it 
can be parametrized as n = (sin#, 0, cosO) which trans- 
forms eqn (0) into a free energy per unit area 



^33 
z .10 



dz 



-W cos 2 O , (4) 



where f 13 (9) = 1 - <5sin 2 6, 8 = {K 33 - Ku)/K 33 , 6 = 
9(z = 0). 

The absence of explicit ^-dependence in the free energy 
(0) implies the first integral 



Boundary conditions read 
09 



= const. 



K 33 fi 3 (9 







1. 



z = L) 



-W sin 20 o , 



(5) 



(6) 



Here we assumed that the director angle at the boundary 
z = L is fixed. 

Integrating eqn (^), together with the boundary con- 
ditions, yields 

E{6, 8) = E(9o,5) + [E(9 L ,S) - E{9 ,S)]z/L, (7) 
[E(9 Ll S)-E(9 0l S)Wf 13 (9o) = (L/2A)sin20 o , 



where E(9, S) = J Q yj f\ 3 {x)dx is the incomplete elliptic 
integral of the second kind, and A = -K33/W is the an- 
choring extrapolation length. For small angles 0q eqn (]?]) 
can be simplified and has the form 



E(0,8)=E{e L ,S) 



z + A 



(8) 



Note, that for small angles 9l and, correspondingly, for 
small 9, E(8,8) = 9 and we have linear dependence of 
the director angle on the z coordinate 



8n x (z) w 9(z) = 



z + A 
L + X' 



(9) 



Using eqn (||), one can fit the director profiles 9(z) with 
A = K 33 /W and 8 = (K 33 — Kn)/K 33 as adjustable pa- 
rameters. To simplify the procedure, it is more appropri- 
ate to fit z{ff) rather than 8(z). 



III. ONSAGER APPROACH 

The Helmholtz free energy in the Onsager approach 
is expressed in terms of the single-particle density, p(l), 
where (1) = (7*1, Jli) represents both position r\ and 
orientation fii. It has the following form [jl3], |l4j 

= J P(l) {lnp(l)A 3 - 1 + 0U(1)} d(l) - 
1 



/(l,2)p(l)p(2)d(l)d(2). 



(10) 



Here j3 = l/fceT, A is the de Broglie wavelength, /i is 
the chemical potential, U is the external potential energy 
(including the surface potential), and /(l, 2) is the Mayer 
/-function 



/(l,2) = exp [-V(l,2)/fc B T]-l, 



(11) 



where elongated particles interact pairwise through the 
potential V(l, 2). 

The equilibrium single-particle density that minimizes 
the free energy (|l0|) is a solution of the following Euler- 
Lagrange equation 

lnp(l)A 3 -/3 A1 + /3C/(l)-y /(l,2)p(2)d(2) = 0, (12) 

which can be obtained from the variation of the func- 
tional (|Io|). In practice, we find it more convenient to 
minimize the functional (|io| ) instead of solving the inte- 
gral equation (12). 



A. Bulk problem 

In the bulk problem, the single particle density is in- 
dependent of position, p(l) = p(Qi). Then the integrals 
over position may be performed directly. To perform 
the integration over the angles, we expand the Mayer /- 
function and the single-particle density in spherical har- 
monics 



/(I, 2) 



/^.^(r ia )*^A- < -(ni > n 3> n t .) 



(13) 
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(14) 



where $ ei ' ta ' tr (Sli, l~i 2 , tl r ) is a rotational invariant |l5|] 

^— ' V mi m<2 m r J 
>£ imi (Oi)^ 2m2 (f2 2 )Q rmr (f2 r ) . (15) 

Here, ( e e> , ) is a 37 symbol, J~i r is the direction of 
the unit vector T\ 2 = r*i2/Vi2 where r%2 = r% — r 2 , Q m 
is a reduced spherical harmonic. In writing eqn (|l4|), we 
assumed that the director n is pointing along the z axis. 

To minimize the grand potential it is convenient to ex- 
pand the logarithm of the density in spherical harmonics 



even 



(16) 



The grand potential then can be rewritten in the form 

n-i \ 2ltV u , ^ 

— = -v47r(l + pn)po + °epe + j= \ = Pepe, (17) 



VW+T' 



where the coefficients pt may be expressed in terms of 
the parameters ci, and where the pair-excluded volume 
is expanded in coefficients 



4,t / r 2 drf m {r). 



(18) 



The grand potential was numerically minimized with re- 
spect to variation of the parameters C£, by the conjugate 



gradient method Jig] . The resulting single-particle den- 
sity was used to calculate elastic constants, for different 
values of the order parameter. 



B. Elastic constants 



To evaluate the elastic constants we used the expres- 
sions derived by Poniewierski and Stecki 17, [l8| Ell 



Kn — M xxxx — M; 



# 



yyyyi 



(19) 



22 



M XX yy Myy X 



yyx.c- 
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33 
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M 
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where 



^k B T J dri2dnidn 2 /(r 12 ,ni,n2) x 

r 12 r 12/ '( COS ^l)/°'( COS ^2)^l7^25- 



(20) 



All integrals are evaluated in a local coordinate frame 
with the z axis parallel to the director at point r. 

As discussed in appendix |a|, performing the integra- 
tions over the angles and making use of the properties of 
3j symbols and spherical harmonics, one obtains a sim- 
plified expression previously given by Stelzer et al 1201] : 



/3#u = #(1,1,3), 
(3K 22 = #(1,-5,-1), 
/3# 33 - #(-2,4,-4), 



(21) 



where 



J 



4tt 



2 even 



K( ai ,a 2 ,a 3 ) = — Y,i(e+i) p: 



CCA) 



6 ai + a 2 £(£+l) 



03 /3 (€ + 3)1 / (21)! 
5 V 2(£-l)\ V (2* + 5)! 



I 

V21+ 1 5 

' „ „ r«+2.2 

I 



(21-2)! ^,2 
(2£ + 3)! 



and where I ll > l2 > 1 are radial integrals over the expansion 
coefficients of the Mayer /-function. 



C. Slab geometry 

In slab geometry, with the z direction normal to the 
surfaces, the single-particle density and the external po- 
tential depend on the z coordinate only. To perform the 
integrations over the angles, we expand the Mayer /- 
function into rotational invariants, similar to the bulk 



system (13). The single-particle density and its loga- 
rithm can be also expanded in spherical harmonics 



p(z, il)=^2 pe m (z)Y e * m (Sl), 
In p(z, ft) = } j ce m (z)Ye m (fl). 



(22) 



The difference from the bulk case, eqns ( |14| , |16[ ), is that 
the director is allowed to vary in the xz plane, so an 
expansion in Legendre polynomials (to = 0) is not suf- 
ficient. Conducting angular and x- and y-integrations 
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gives 



-J dzV iir(l + (3 fi) poo (z) + 

Y^e,m( c t™{z) + f3Ui m (z))p£ m + 
2ir J dzidz 2 Ai u i 2im (zi2)p£ im (zi)pi 2m (z 2 ) 



(23) 



where the pair-excluded area at given z-separation 212 is 
expanded in coefficients 



i x t 2 1 

m m 



A£ 1: l 2: m{zi2) = -27T ^ 

sds ft u t 3t t(r 12) Pi(cos 





(24) 



Here s 2 = (xi - x 2 ) 2 + (yi - y 2 ) 2 , tan6» r = s/(zi - z 2 ), 
rn = —m. Integrating, we took into account that the 3j 
symbol vanishes unless mi + m% + m r — 0. 

To obtain equilibrium single-particle density profiles, 
the grand potential then was tabulated on a regular grid 
of points Zi and numerically minimized with respect to 
variation of the parameters Q m (zi), by the conjugate gra- 
dient method. 



D. Anchoring energy 

To obtain the microscopic expression for the extrapola- 
tion length A, we start from the equation for the single- 
particle density Jl^). Assume that the solution for a 
ground state, i.e. for homeotropically aligned liquid crys- 
tal in slab geometry, is given by the single-particle den- 
sity po- Consider a small perturbation around the ground 
state, p = po + Sp. To first order in dp, eqn ( |l2| ) can be 
written as 



Po(l) 



/(l,2)5p(2)d(2) 



(25) 



In the case of slowly- varying director fields we assume 
that the free energy functional is locally in equilibrium. 
This is equivalent to the mathematical simplification 



p(r,ft) = po (r,n(r) ■ O) . 



(26) 



Then the density perturbation can be written in terms 
of the perturbation of the director 



6p(r, ft) = p' (r, no ■ ft) 5n(r) ■ ft, 



(27) 



where the prime denotes a partial derivative with respect 
to (n ' ft)- 

In slab geometry, with the z axis normal to the sur- 
faces, the single particle density po depends on the z 
coordinate only and can be expanded in spherical har- 
monics similar to eqn ((2^). Note, that po does not de- 
pend on (j> which implies m = in expressions (53). 
We also assume that the director is parallel to the xz 
plane, 5n = (5n Xl 0,0). Conducting angular integrations 
in eqn (|2^) and making use of the properties of 3j sym- 
bols we obtain 



Cf 



1 (zi)Sn x (zi) = 47T 



U 2 (e 2 



2=2 



Ai u i 2t i(z2 - z 1 )pe 2 (z 2 )6n x (z 2 )dz 2 . 



(28) 



Equation (|28|) is a homogeneous Fredholm equation of 
the second kind. It allows one to calculate the direc- 
tor profile (for small director deviations from the ground 
state, no = e z ) once the single-particle density po of the 
ground state is known. Eqn ( p^ ) is not valid for every £1, 
in spite of the derivation. This is because we are trying 
to map the single-particle density variation onto the di- 
rector variation (eqn |26|). However, this equation should 
be valid for the leading term in the density expansion, 
£i=2, which we consider below. 



First, we construct the solution to the eqn (|2§|) in the 
cell bulk. The kernel Ag 1 ^ 2 ^{z 2 — z\) is a short-ranged 
function. It equals zero for \z 2 — Z\\ > a, where a is the 
molecular length. The bulk director is a slowly-varying 
function on this length scale. Therefore we can expand 
it in a Taylor series 



8n x (z 2 ) = 8n x (zi) + 8n' x {zx)zx 2 + 5n"(zi)z 2 2 /2 + 



(29) 



and restrict ourselves to second order in the expansion. 
Then the equation for the director (^|) can be- rewritten 
as a second-order linear differential equation 



a 2 (z)8n'^.(z) + ai(z)Sn' x (z) + ao(z)8n x (z) = 0, (30) 



where 



a n (z) 



A 2: i 2: i{z 2 - z)pt 2 (z 2 )^dz 2 ~ c 2 (z)8 n , . 



(31) 
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In the bulk, the expansion coefficients c 2 and pe 2 do not 
depend on the z coordinate. Expansion coefficients of the 
excluded area, Ai li £ 2i £ r (z), are even functions of z, which 
immediately implies that a\{z) — in the bulk, since the 
integrated function zAi lt i 2 j r (z) is odd (because a\{z) is 
zero in the bulk, we expanded Sn x (z) up to second order). 



J 



To prove that ao(z) — in the bulk, we consider again 
the equation for the single particle density (|l2|). Perform- 
ing the usual expansion of density and /-function (13 I4| ) 
and integrating over the angles and x,y coordinates we 
obtain 



ca(zi) 
cejzi) 



4x011- 4tt A ,oA z 2- zi)po(z 2 )dz 2 , 



even „£, 

-4tt ^2 / A ei j 2 . (z 2 - zi)p e2 (z 2 )dz 2 , (h^O) 

I 



(32) 



where Ai 1 j 2 _ x are the expansion coefficients of the ex- 
cluded area in a series of spherical harmonics. It is 
easy to show that the left hand side of the second equa- 
tion in (J3|) equals a (z) in the cell bulk. Indeed, us- 
ing the symmetry of the excluded volume expansion 
coefficients we can write — 4n J Q L A 2 j 2 ,i(z 2 — Zh)dz 2 = 
47r J^° oo A 2 ^ 2 fi(z)dzS 2 j 2 = V 22 S 2y £ 2 , which converts the 
expression for ciq(z) to the left hand side of the second 
equation in (|32"|). In fact, the conclusion ciq{z) = in the 
cell bulk is a consequence of the invariance of the grand 
potential with respect to director rotations. 



Therefore, eqn ( p0[ ) in the cell bulk reduces to 5n"(z) = 
which means we have a linear dependence of the di- 
rector on z coordinate, in agreement with the result of 
elastic theory @. This also means that, in the cell bulk, 
cq + c\z is an eigenvector of the Fredholm equation (28). 



To solve eqn (|28|) near the cell boundary is a much 
more challenging task. Numerically it can be done using, 
for example, an iterative method ]l6| |. Here we try to 
construct a crude analytical solution which will give us 
some qualitative understanding of what is happening in 
the interface region. 



To begin with, we simplify eqn ( 28 ) . The sum over £ 2 in 
the right-hand side of eqn ( p8f ) is converging very fast (for 
ellipsoids with the elongation e = 15, every term is about 
ten times smaller than the previous one). Therefore, we 
truncate the sum leaving only the t 2 = 2 term. Then, 
the kernel of the Fredholm equation A 2%2 _\(z\ — z 2 ) can 
be expanded in a Taylor series 



A 2 . 2 , 1 {z 1 - z 2 ) = J2 Z "^ 

n=0 



(*2) 



(33) 



where $ n (z) are some functions. A practical example 
of such an expansion is given in appendix |b|. The ker- 
nel is then separable, and the problem is reduced to the 
solution of a set of linear algebraic equations. Indeed, 



substituting eqn ( |33D into eqn (|28|) we obtain 

oo 

5n x (z) = C2 1 (z)J2b n z n , (34) 

where the b n are some constants. Substituting (|34|) back 
into eqn ( p8| ) we obtain an infinite set of linear equations 
for the coefficients b n . To obtain an analytical expression, 
we perform further simplifications. First we note that, in 
the cell bulk, the director is a linear function of the z 
coordinate. Therefore, to a good approximation, we may 
retain only the first two terms in eqn (^8|) since in the 
cell bulk c 2 (z) = const. Using again the director profile 
given by elastic theory, eqn (|9|), we obtain an expression 
for the extrapolation length 



x= bo = 4tt Jq 00 zA 2:2:1 (z)p 2 (z)c 2 1 (z)dz 
h 1-4tt J °° A 2 ^ 1 (z)p 2 {z) C 2 1 (z)dz' 



(35) 



This expression is able to give qualitative explanations 
of the anchoring phenomenon in our system. In the 
ideal case, often considered in phenomenological ap- 
proaches p9| , the density and order parameter are as- 
sumed to be constant in the cell, i.e. p 2 {z)c^ x {z) = 
const. The anchoring appears only due to the pres- 
ence of the interface, which breaks translational sym- 
metry. According to eqn (|32"|), c 2 = —V 22 p 2 , therefore 
47T Jq 00 A 2i2 ^{z)p 2 {z)c2 1 dz = 1/2. Then the anchoring 
coefficient is proportional to the first moment of the ex- 
cluded area coefficient A 2 ^ 2 _i(z) 



A 



8tt 



zA 2j2t i(z)dz, 



(36) 



and does not depend on the value of density or order pa- 
rameter, since excluded area, as well as excluded volume, 
are completely defined by the geometry of the overlap- 
ping molecules. This means that in the ideal case of a 
uniform ncmatic, the anchoring coefficient W = K^/X 
has the same dependence on the order parameter as the 
elastic constant K 33. 
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In reality, one observes rather strong subsurface vari- 
ations of the density and order parameter. The ratio 
Pi{z)lc-i{z) then comes into play and contributes to the 
overall dependence of the extrapolation length on the or- 
der parameter. The numerical r esults which show this 
dependence are presented in sec. VI B. A simple physi- 



cal explanation is also possible. Presmectic ordering and 
higher density of the nematic phase in the interface re- 
gion indicate that the nematic liquid crystal is at different 
state point. The mesophase is more ordered at this state 
point and this ordering is less sensitive to the density 
variation in the cell bulk. Since this ordering defines the 
director profile at the interface, it also affects the depen- 
dence of the extrapolation length on the state point. 



IV. THERMAL FLUCTUATIONS 

Another simulation method to measure bulk elastic 
constants Kn and the zenithal anchoring energy W is 
based on the measurement of the director fluctuation 

EaL M E3. 



amplitudes Sn in the liquid crystal cell 22 



Consider again slab geometry with homeotropic anchor- 
ing of the director at both cell surfaces. Consider a small 
perturbation of the director around the equilibrium dis- 
tribution 



n(r) = e z + Sn(r). 



(37) 



Minimizing the free energy (|2[) and linearizing the equa- 
tions for the director and boundary conditions with re- 
spect to Sn, we obtain 



Ku(5n X:XX + 8n Vtyx ) + K 2 i{Sn x ^ y - 5n Viyx ) + K 33 5n X:ZZ = 0, 
Kii(5n XtXy + Sny^yy) + K 2 2{5n x ^ x - Sn y , yx ) + K 33 Sn x , zz = 0, 



d 

WSn + K 33 —Sn 

oz 



= 0, 



z=L 



d 

WSn — A'33— Sn 

Oz 



= 0. 



z=0 



(38) 



Here we adopt the notation Sn a ,pi — d/3dj(5n a ). 

We now expand Sn(r) in a series of orthogonal func- 
tions 



Sn(r) = i ]T e 1 ' 

q±,q~ 



Sn^\q^q z )e^ + Sn^ (q ± , q z )e~^ 



(39) 



density 



Q a p{r) 



V 

N 

3 



i 

1 „ 



!ap = 2 I u ia u il3 ~ ^af} 



where q± = (q x ,q v ), and <5n ( ) = (ix - £)/ (ix + 
to satisfy the boundary conditions. Here we have in- 
troduced the dimensionless wave vector \ = QzL and 
anchoring parameter 



WL 

w 



L 
A ' 



(40) 



where A is the extrapolation length Q . The wave vectors 
q z form a discrete spectrum because of confinement in 
the z direction which depends on the anchoring energy 
W . The explicit form of the q z spectrum is given by the 
secular equation: 



(£ 2 -X 2 )sin X + 2£xcosx = 0. 



(41) 



In molecular simulations, rather than measuring direc- 
tor fluctuations, it is more convenient to measure fluctu- 
ations of the second-rank order tensor components (fol- 
lowing Ref. p6[). We define the real-space order tensor 



where a,/3 = x, y, z, in terms of the orientation vectors Ui 
of each molecule i (we consider only uniaxial molecules). 
If we assume that there is no variation in the degree of 
ordering, and neglect biaxiality of the order tensor, we 
may write 



3 1 
)afi(r) = -Qn a (r)n (r) - -Q5 af} , 



where Q is the order parameter, i.e. the largest eigenvalue 
of Q a p(r). 

Measurements are performed directly in reciprocal 
space. The Fourier transform of the real-space order ten- 



Q a0 (k) = J Q a p(ry kr dr = Q^e*" 
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Then the fluctuations (\Q a p{k)\ ) can be easily mea- 
sured from simulations 



\Qafi(kz)\ = 



(42) 



and compared with the prediction of elastic theory 

,2 i <*2 



X 



K 33 ^ ql{2i + X 2 +e\ 



X 



K + X 



ix - 



-A c'^-x) _ i 



K-X 



(43) 



where (...) denotes an ensemble average, k = k z L. 

We measure Q and (iQazikz)^ from simulations, 

eqn ([42]), and then compare with the theoretical pre- 
diction, eqn (|43|), which is parametrized by L, A and 
-K33. Both the permitted q z spectrum, and the variation 



of \ \Qaz(k z )\ j with k z , are sensitive to the anchoring 
strength parameter £ = L/X . 



MOLECULAR MODEL AND SIMULATION 
METHODS 



We performed Monte Carlo (MC) simulation of a liq- 
uid crystal confined between parallel walls, with finite 
homeotropic anchoring at the walls. We used a molecu- 
lar model which has been studied earlier in this geometry 
p5| . The molecules in this study were modeled as hard 
ellipsoids of revolution of elongation e = a/b = 15, where 
a is the length of the major axis and b the length of the 
minor axis. Such systems cannot form smectic or colum- 
nar phases |27|, and the phase transitions are not ther- 
mally driven as they are for most mesogens. Therefore, it 
should be borne in mind that simulation results are model 
specific, however the advantage of using systems of hard 
ellipsoids is that they exhibit a simple phase behaviour 
with some resemblance of that of real nematogens. In ad- 
dition, as the elongation increases, the nematic-isotropic 
phase boundary approaches the Onsager limit. This 
eases the comparison between simulation and density- 
functional Onsager-type theories. 

The phase diagram and properties of this family of 
models are well studied §7], |?§ ||, |?| It is useful 
to express the density as a fraction of the close-packed 
density p cp of perfectly aligned hard ellipsoids, assuming 
an affinely-transformed face-centered cubic or hexagonal 
close packed lattice. For this model, temperature is not 
a significant thermodynamic quantity, so it is possible to 
choose k-gT = 1 throughout. 



The slab geometry is defined by two hard parallel con- 
fining walls, which cannot be penetrated by the cen- 
ters of the ellipsoidal molecules. Packing considerations 
generate homeotropic ordering at the surface, since the 
molecules have more free volume if their axes are normal 
to the wall. Surface anchoring in this system has been 
studied recently by applying an orienting perturbation 
at one of the walls and observing the response at the 
other [ p5| and by measuring amplitudes of the director 
fluctuations |^2| . This yielded an estimate of the extrap- 
olation length A « 17.56 « 1.16a at one particular state 
point corresponding to the value of the order parameter 
Q w 0.85. 

A sequence of runs was carried out for systems of 
N = 2000 particles using the constant-TWT ensemble, 
allowing typically 10 5 MC sweeps for equilibration and 
4 x 10 7 sweeps for accumulation of averages (one sweep 
is one attempted move per particle). The wall separa- 
tion was fixed and equal to L z w 4.93a. Note, that the 
size of the simulation box L z is not equal to the liquid 
crystal cell thickness L appearing in the elastic theory. 
The difference can be ascribed to partial penetration of 
the walls by the liquid crystal molecules, and/or forma- 
tion of a solid layer near the walls. It is possible to write 



L = L, 



2L W , where L w characterizes the wall. In our 



previous publication we found L w ss 0.3a p2| . 

The same molecular model and interaction of the 
molecules with the walls was adopted for the Onsager cal- 
culations. It was found sufficient to include terms with 
< / < 8 in the expansion of In p, while terms with 
< I < 10 were included in the expansion of the pair 
excluded volume coefficients Vu and the excluded area 
coefficients Ao, o n m . 



VI. RESULTS AND DISCUSSION 



A. Onsager theory, bulk 



Minimization of the free energy functional (|10|) was 
carried out for several values of chemical potential \i. 
From the single-particle density we evaluated Frank elas- 
tic constants, which, together with the values of the bulk 
density as a fraction of the close-packed density and the 
value of the order parameter are given in table Q. The 
results are typical for the elastic constants of prolate bod- 
ies. They increase with fluid density (order parameter) 
and K Z >K 1 > K 2 |TJ. 

It is essential to carry out bulk calculations if we want 
to know both anchoring extrapolation length A = ^33/!^ 
and anchoring strength W . The bulk problem provides 
us with the absolute values of the elastic constants; the 
elastic theory includes only ratios of elastic constants in 
the expression for the director profile. 
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TABLE I: Onsager calculations in the bulk. Reduced den- 
sities pj pep, values of the nematic order parameter Q, and 
elastic constants for hard ellipsoids with elongation e = 15. 





pl Pep 


Q 




K22 


#33 


1.2 


0.0293 


0.7193 


0.7951 


0.2925 


3.5290 


1.3 


0.0303 


0.7532 


0.8883 


0.3294 


4.4024 


1.4 


0.0313 


0.7779 


0.9677 


0.3614 


5.2387 


1.5 


0.0323 


0.7973 


1.0398 


0.3909 


6.0704 


1.6 


0.0332 


0.8132 


1.1076 


0.4190 


6.9113 


1.7 


0.0341 


0.8266 


1.1726 


0.4463 


7.7694 


1.8 


0.0349 


0.8382 


1.2361 


0.4731 


8.6496 


1.9 


0.0358 


0.8482 


1.2985 


0.4997 


9.5547 


2.0 


0.0367 


0.8571 


1.3605 


0.5263 


10.4883 


2.1 


0.0376 


0.8650 


1.4225 


0.5531 


11.4510 


2.2 


0.0384 


0.8721 


1.4848 


0.5801 


12.4448 


2.3 


0.0393 


0.8785 


1.5477 


0.6076 


13.4716 


2.4 


0.0402 


0.8843 


1.6113 


0.6354 


14.5318 


2.5 


0.0410 


0.8896 


1.6758 


0.6639 


15.6268 



0.90 
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FIG. 1: Order parameter Q as a function of the density (mea- 
sured far from the surface). The density is expressed relative 
to the closed-packed density p cp for perfectly aligned ellip- 
soids. Filled circles: simulation results. Open circles: On- 
sager theory calculations in slab geometry. Solid line: On- 
sager theory, bulk. 



B. Onsager theory, slab geometry 

Minimization of the grand potential was carried out in 
slab geometry for the same values of chemical potential 
fi as considered in bulk. From the density and order 
parameter profiles we were able to extract values of these 
quantities in the central part of the cell which agreed with 
the bulk Onsager calculations. 

Together with the simulation results, we plot depen- 
dences of the bulk order parameter versus density in 
Fig. |l|. The Onsager theory does not describe the bulk 
equation for the state perfectly: the predicted bulk den- 
sity for a given order parameter value is larger than the 
value obtained in the simulation. 

With this choice, the same values of the chemical po- 
tential and the same slab dimensions, minimization of 




20 40 60 80 100 120 140 



z/b 

FIG. 2: Typical profiles of director angle 9(z) in the slab 
geometry. Orienting fields are applied in the region 100 < 
z/b < 148 near the right wall, favouring director angles of 
7r/4 and 7r/2 relative to the surface normal. The left wall is 
unperturbed. Solids lines: Onsager theory. Circles: results of 
fitting the profiles in the bulk region with the prediction of 
elastic theory. 
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FIG. 3: Ratio of the elastic constants, S = (-K33 — Ku) / K-^. 
Circles: Poniewierski - Stecki expressions. Squares: Onsager 
theory in slab geometry, with wall anchoring field at angle 
a = 7r/4, fitted with the results of elastic theory. Triangles: 
Onsager theory with a = 7r/2, fitted with elastic theory. 



the grand potential was carried out for the system with 
an external field applied near the right wall. Profiles 
of director angle 0(z) are compared with elastic theory 
(eqn (||)) in Fig. 0. The elastic theory has been fitted 
to the director angle profiles predicted by the Onsager 
theory using two adjustable parameters, the extrapola- 
tion length A and the elastic constant ratio 5. Note that 
only part of the bulk region has been used for fitting, 
20 < z/b < 80 where the elastic theory is applicable. 

The anisotropy of elastic constants S obtained from 
fitting is shown in Fig. |j, together with the results of cal- 
culations using the Poniewierski and Stecki expressions. 
The value of 5, and hence the splay constant K\, comes 
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into play only as the deformation angle 6{z) increases. 
Therefore, for the external field with easy axis 6l = 7r/4, 
the error in the determination of S from simulation data 
is quite large. 

The dependence of the extrapolation length, A, on the 
order parameter Q is shown in Fig ^. It is seen that 
Onsager theory predicts the extrapolation length to grow 
linearly with the order parameter, which is completely 
different from the experimental results in thermotropic 
liquid crystals, where A decreases with the increase of 
the order parameter as Q~ 2 . 

We have also carried out minimization of the grand 
potential for the system without an external field. As 
a result, we obtained the single-particle density with a 
homeotropic distribution of the director in the cell. The 
dependence of the extrapolation length on the order pa- 
rameter was then calculated using eqn (|3^) and is also 
shown in Fig [|. The results qualitatively agree with the 
results of fitting obtained by combining Onsager theory 
and elastic theory. The extrapolation length tends to 
grow with increase in the order parameter and has the 
same order of magnitude. Plugging this single-particle 
density into all equations preceding eqn (|3^) we were 
able to check the approximations we did deriving this 
equation. We found that the most crude approximation 
is the truncation of the sum in eqn (|34|). This is not 
a fundamental problem and can be easily corrected by 
taking into account a sufficient number of expansion co- 
efficients. However, this points out that the dependence 
of the director on the z coordinate near the cell surface 
is different from the linear dependence in the cell bulk. 
Moreover, to obtain correct quantitative values of the an- 
choring coefficient, or extrapolation length, one needs to 
know the director distribution at the surface. The as- 
sumption 5n = const in the interface region, which has 
been made in ||, may lead to absolutely incorrect esti- 
mates of the anchoring coefficient. 



C. Simulation 

Simulations were carried out in slab geometry for sev- 
eral values of the number density. The density variation 
0.25 < p/pc P < 0.34 provides a sufficient range of order 
parameter variation in the nematic phase, 0.68 < Q < 
0.88 for us to study the effect of ordering on elastic be- 
haviour. 

The order tensor fluctuations in reciprocal space were 
calculated using expressio n j42| ). To fit the simulation 
results we used expression Q43Q. Recall that the size of the 
simulation box L z is not necessarily equal to the liquid 
crystal cell thickness L appearing in the elastic theory. 
We found that L — L z + 2L W , with L w sa 4.56 = 0.3a, 
almost independent of the density. 

Using this value of L w we adjusted the extrapolation 
length A and the elastic constant K33 to obtain the best 
fit to the fluctuation amplitudes. Typical order fluctua- 
tion amplitudes together with the corresponding fitting 




0.7 0.75 0.8 0.85 0.9 

Order parameter, Q 



FIG. 4: Extrapolation length A as a function of the order 
parameter. Circles: Onsager theory in slab geometry, direc- 
tor profiles fitted with elastic theory, for anchoring field with 
a = 7r/4. Triangles: the same, but for a = tt/2. Diamonds: 
Onsager theory in slab geometry, with no field, extrapola- 
tion length calculated using equation (|35|). Filled squares: 
Monte Carlo results obtained by measuring director fluctua- 
tions. Open squares: Monte Carlo results, with applied field 
near the right-hand wall. 




FIG. 5: Fluctuations of the director (arbitrary units) as a 
function of wave-vector (normalized by the molecular minor 
axis length b). Symbols: Monte Carlo results. Solid lines: 
elastic theory, fitted to parameters discussed in the text. In- 
set: fluctuations multiplied by (k z b) 2 to emphasize structure 
at higher wave numbers. 



curves are plotted in Fig From this plot one can see 
that the fitting curves agree quite well with the simula- 
tion results for small values of the wave- vector k z . At 
higher k z , the structure is not perfectly reproduced, as 
one would expect for a theory valid for long-wavelength 
fluctuations, but the agreement is satisfactory. 

The elastic constant, K^b/k^T versus order parame- 
ter Q is plotted in Fig [| The agreement with the Onsager 
theory is satisfactory, especially if we take into account 
that the equation of state is not perfectly reproduced by 
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FIG. 6: K^b/ksT. Circles: Onsager bulk calculations. 
Squares: MC simulation results. 



the Onsager theory. 

The dependence of the extrapolation length, A, on the 
order parameter Q is shown in Fig [|, together with the 
results from the Onsager theory. It should be noted that 
combining the elastic approach with the Onsager calcu- 
lations does not allow one to determine, separately, L w 
and A. Therefore, the results of the Onsager theory in 
Fig U really represents A + L w , which is one of the possi- 
ble origins of the systematic difference between the two 
approaches. 

To double check the results obtained by examining 
the director fluctuation amplitudes, we performed the 
same type of experiments as in the Onsager slab sys- 
tem. Within a range 7.56 of the right-hand wall, a 
strong coupling field was applied to molecular orienta- 
tions, V cxt ~ (ui ■ e 2 ) 2 , aligning the molecules near the 
right wall parallel to it. After the system was equili- 
brated, the director profile was fitted to the result of the 
elastic theory, eqn (||). The dependence of the extrapola- 
tion length A on the order parameter Q is also shown in 
Fig " 
met 



One can see that the agreement between the two 
lods is quite good. 
Finally, we plot the dependence of the anchoring en- 
ergy coefficient, W — -K33/A, which is an intrinsic 
charachteristic of the interface region, in Fig. 0. For On- 
sager theory, X(Q) is given by the results in slab geom- 
etry fitted with the elastic theory (Fig. ||), K 33 (Q) - by 
the Poniewierski - Stecki expressions (Fig. ^|). For MC 
simulations, we used the elastic constant obtained from 
the analysis of the director fluctuation amplitudes. All 
methods predict that the anchoring coefficient is a non- 
monotonic function of the order parameter, even though 
the actual variation is small. 



VII. CONCLUSIONS 

We have studied the dependence of the zenithal sur- 
face anchoring coefficient on the order parameter of a ly- 
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FIG. 7: Anchoring energy W as a function of the order pa- 
rameter. Circles: Onsager theory (with a smooth curve to 
guide the eye). Open squares: MC results obtained by mea- 
suring director fluctuations. Filled squares: MC results, field 
near the right-hand wall. 



otropic nematic liquid crystal modeled by hard ellipsoids. 
Several techniques have been used: Onsager theory com- 
bined with elastic theory; Monte Carlo simulations fitted 
to elastic theory; analysis of the director fluctuation am- 
plitudes obtained from Monte Carlo simulations. The 
results of these methods agreed qualitatively with each 
other. 

Perhaps the most interesting aspect of this study is the 
increase of the anchoring extrapolation length with the 
increase of the nematic order parameter. This implies 
that the bulk elastic moduli in our system grow faster 
than the surface anchoring strength. This is opposite to 
the experimentally observed behaviour in thermotropic 
nematics, where the extrapolation length goes down with 
increase of the order parameter. 

A microscopic semi-qualitative expression for the ex- 
trapolation length allowed us to conclude that subsurface 
variations of the single-particle density, mainly defined 
by the nematic order parameter and density variation, 
contribute substantially to the anchoring phenomenon. 
We showed that for an ideal system, in which the single- 
particle density in the cell is assumed to be uniform, the 
extrapolation length does not depend on the nematic or- 
der parameter. This dependence is therefore associated 
with the subsurface variations of the single-particle den- 
sity. 

We now turn to a brief discussion about possible gen- 
eralizations of this work. First, it would be interesting to 
go beyond the limits of Onsager theory and use the direct 
pair correlation function in the nematic liquid inst ead o f 
the Mayer /-function. Second, as was shown in sec. [II D, 



the anchoring coefficient can be calculated if we know 
the single-particle density and the direct pair correlation 
function (or Mayer /-function in case of Onsager theory) . 
This has also been done by Stelzer et al j^] with numer- 
ous approximations. Combining elastic theory and local 
density functional theory one can avoid these approxi- 
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mations, or at least check their validity. This work is in 
progress. 

Acknowledgments 

This research was supported by EPSRC grants 
GR/L89990 and GR/M16023. D.A. acknowledges the 
support of the Overseas Research Students Grant; 
M.P.A. is grateful to the Alexander von Humboldt foun- 
dation, the British Council ARC programme, and the 
Leverhulmc Trust. 



Appendix A: Elastic constants 

Performing the integrations over the angles in expres- 
sion (|2(i| ) and making use of the properties of 3j symbols 
and spherical harmonics one obtains 

PK U = «(1, 1,1, 1,1,1), (Al) 
[3K 22 = k(1, -1,1, -1,1,-1), 
/3K 33 = «(-2, 0,-2, 0,-2,0), 

where 



K(b!,b2,b 3 ,b4,b 5 ,b 6 ) 





1)! 


h 


V- 




5 


(i+ 


3)! 


'b 5 


V- 




5 



1 



I 
-1 



£ £-2 2 
1 -1 

t £ + 2 2 
1 -1 



b 4 /3 



5 



1 -1 1 



^ ^-2 2 
1 1 -2 

£ £ + 2 2 
1 1 -2 



h. 

5 



'J+2,2 



£ £ 2 
1 1 -2 



+ 



(A2) 



and J*i>4>>* are radial integrals over the expansion coeffi- 
cients of the Mayer /-function 



(A3) 



Expression (Al) is the same as obtained before by Stelzer 
et al (20), except that we used 3j symbols instead of 
Clebsch-Gordan coefficients and a different normalization 



of the single-particle density. To simplify (A2) we use the 
relation between the 3j symbols 



3f £ £ + 2 2 
2 I 1 1 -2 



1(1 £+2 2 
2 I 1 -1 



(A4) 



and take into account that sums with 63, 65 and 64, bg 
are identical, since 



li I2 2\_ ,_ r ,i 1+ e 2 ( £2 £1 2 \_ (£ 2 £1 
1 -1 ) ~ [ > \ -1 1 0) ^1-1 



2 

-1 



and 7^1^2,2 _ j« 2 ,*i,2_ Performing simplifications and 
taking into account explicit expressions for the 3j coeffi- 
cients we obtain expression (|2l|). 



Appendix B: Expansion of the kernel 

From the numerical data we found that the kernel of 
the integral equation A221 (z) can be accurately approxi- 



mated with a Gaussian function 

1^22 



12,2,1 



exp -(W 



(Bl) 



where we took into account that — 4tt J_ A2.2, \{z)dz — 
V%%. Here I is a geometrical parameter which depends 
only on the elongation of the molecules. 

Using the generating function for the Hermite polyno- 
mials [fc>3fl 



00 t n 
exp(-t 2 + 2tx) = V H n {x) — 



(B2) 



n=0 



one can write the kernel in a separable form 

V22 



42,2,1(22 - zi) 



4^ GXPh(Z2/0JX 

OO n 



(B3) 



n=0 



similar to the general case we used in sec. JIID. 

Using this approximation to the kernel one can show 
that the expression for the anchoring coefficient (|3^ ) 
reads 



(B4) 



For ellipsoids with elongation e = 15, we found I « 5.466. 
This results in an anchoring coefficient A « 3.16 which is 
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much lower than the actual value of the anchoring in the 
system with the density and order parameter variations 
at the surfaces. Hence, changes in the density, order 



parameter, and, as a result, in the director profile in th e 
interface region cannot be neglected (see also sec. VT Bp . 
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